Theoretical study of d-wave superconductivity in doped bilayer graphene 
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We introduce a microscopic model on the honeycomb bilayer, which in the small-momentum limit 
captures the usual description of bilayer graphene. Due to assumed strong interlayer hopping it 
reduces to an effective honeycomb monolayer model with also third neighbor hopping. We study 
interaction effects in this effective model focusing on possible superconducting instabilities. We 
find d x 2_ y 2 superconductivity in the strong coupling limit of an effective t J- model-like description 
that gradually transforms into d + id time reversal symmetry breaking superconductivity at weak 
couplings. In this limit the small momentum order parameter expansion is (k x + ik y ) 2 (or (k x — 
iky) 2 ) at both points (valleys) of effective, low-energy description. The relevance of our model and 
investigation for the physics of bilayer graphene is also discussed. 



I. INTRODUCTION 

Interaction effects are expected to be important for 
the physics of bilayer graphene and ma y cause a forma- 
tion of correlated many-body phasesP^ This needs to be 
contrasted to intrinsic monolayer graphene in which a 
vanishing density of states at the Dirac points suppresses 
the influence of electronic correlations Recent experi- 
ments on suspended bilayer graphenepH^ which is free of 
substrate effects, reveal a gapped state at and around the 
charge neutrality point. These states may be topolog icaP 
due to the observed^ conductance of the order of e 2 /h 
and may exhibit an anomalous quantum Hall effect, i.e. 
a quantum Hall effect at zero magnetic field. In the most 
recent experiment on high mobility samples from Ref. O 
a completely insulating behavior was found. 

^From the theory point of view, several proposals were 
giverPEDfor the existence of gapped (and gapless) phases 
at the charge neutrality point, including those that break 
the time reversal symmetry. Most of them are based on 
the particle-hole (excitonic) binding which is the most 
natural assumption in the understanding of a gapped 
phase at the charge neutrality point. However, these the- 
ories assume a quadratic dispersion of the electrons in the 
low energy effective description, while direct hopping be- 
tween two sublattices in different layers should lead to 
the linear dispersion ( "triagonal warping" ) that becomes 
relevant at energies < 10 meV. Only if we work at a finite 
chemical potential away from the charge neutrality point, 
that is above ~ 10 meV and below ~ 400 meV where a 
second band becomes relevant,^ the effective quadratic 
dispersion is justified. 

To explore additional possibilities for gapped phases 
in the presence of a finite chemical potential, we discuss 
here superconducting instabilities, especially with an eye 
on the possibility of topological (fully gapped) supercon- 
ductivity on the honeycomb bilayer. This study is moti- 
vated by a phenomenology similar to topological insula- 
tors and topological superconductors, where all transport 
is done via the edge states. The bilayer graphene may be 



also viewed as a strongly-correlated system with a pos- 
sibility to support a layered antiferromagnetic statej^^ 
similar to the Mott physics of high T c superconductors. 
The existence of a layered antiferromagnetic state is sup- 
ported by the most recent experiment with high quality 
samples,- which feature completely insulating behavior 
at the charge neutrality point. 

There is, so far, no systematic study of superconduct- 
ing instabilities in the presence of electron-electron and 
electron-phonon interactions on the honeycomb bilayer 
at finite doping (see, however, Ref. [3T] for fermions in 
the presence of weak electron-electron interactions only 
at zero chemical potential) . To address this question, we 
study in the present paper a microscopic model of a sin- 
gle effective honeycomb monolayer with reduced nearest 
neighbor hopping and third-neighbor hopping, in addi- 
tion to inter-site attractive interactions. The kinetic term 
of the effective model is obtained by integrating out the 
"high-energy" degrees of freedom from the direct inter- 
layer hopping, and the inter-site superexchange interac- 
tion originates from the Hubbard on-site repulsion. This 
model is to a certain degree biased to antiferromagnetism 
and d-wave superconductivity, but preserves the usual 
low-energy description of the bilayer graphene^ and ac- 
counts moreover for the lattice symmetry of the original 
model (the honeycomb bilayer) that may be relevant for 
the symmetry of the superconducting order parameters. 

Our primary interest is to find the most probable sym- 
metry of a superconducting instability on the honey- 
comb bilayer. Furthermore, we aim at an understand- 
ing of the change in the superconducting order param- 
eter and correlations as we go from a monolayer to a 
few-layer honeycomb lattice. This may be of relevance 
for experiments that point to superconducting correla- 
tions in graphite and graphite compounds (see Ref. [2H 
and references therein). The mean-field solution of the 
introduced model yields a time-reversal symmetry break- 
ing d + id-wave superconducting state at weak coupling, 
which quickly but continuously transforms into d x 2_ y 2- 
wave with increasing interaction. Near 3/8 and 5/8 fill- 
ing of the 7r-bands, i.e. near the van-Hove singularity in 



the density of states, the Cooper pairing becomes much 
stronger. Our conclusion is that the d + id superconduct- 
ing instability may be present in the bilayer graphenc 
at finite doping. However, due to the presumed small- 
ness of coupling constant and order parameter, as well as 
strong quantum fluctuations in two dimensions, it may 
be difficult to detect this order experimentally in today's 
graphene samples. 

The remaining part of the paper is organized as follows. 
In Sec. II wc define our effective two-band model on an 
effective honeycomb lattice with third-nearest-neighbor 
hopping and discuss its validity. The model is then, in 
Sec. Ill, solved by a Bogoliubov - de Gennes (BdG) 
transformation for a singlet bond-pairing order param- 
eter, and we discuss the relevant symmetries. Section IV 
presents the phase diagram obtained from a numerical so- 
lution of the BdG equations. In Sec. V, the relevance for 
the physics of the bilayer graphene is discussed, and our 
main conclusions are presented in Sec. VI. Two Appen- 
dices summarize analytically obtained solutions in the 
weak-coupling BCS limit. 



II. MODEL 

The honeycomb bilayer lattice consists of two Bernal- 
stacked honeycomb lattices, each consisting of two tri- 
angular sublattices as illustrated in Fig. 1 such that the 
unit cell contains four lattice sites. The Hamiltonian of 
free electrons on such a lattice is given by 
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Here, the index i — 1, 2 denotes the layer and j enumer- 
ates primitive cells. The sum runs over u — u ,ui,u 2 , 
where u\ = a(§, ^) and H2 = a(§, — ^) are the primi- 
tive vectors of the lattice, and uq=(0, 0) is an auxiliary 
vector for denoting the hopping between sites in the same 
primitive cell. The norm of these vectors is |u| = V3(*, in 
terms of the distance a — 0.142 nm between neighboring 
carbon atoms in each layer, and t ~ 3 eV is the associated 
hopping energy, whereas t± ~ 0.4 eV denotes the inter- 
layer hopping energy, between A sites in two different 
layers. Although the hopping energy between adjacent 
sites B\ and B^ in different layers is estimated^ as ~ 0.3 
eV and thus on the same order of magnitude as t±, it 
yields a fine structure of the dispersion relation (in the 
form of trigonal warping and additional Dirac points with 
linear dispersion) that is only relevant at very low ener- 
gies < 10 meVp2 In the present paper, we neglect this 
fine structure by setting the hopping between the Bl and 
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FIG. 1. (Color online) (a) A view of Bernal stacked honey- 
comb lattices 1 and 2 with corresponding sublattice sites Al, 
Bl, and A2, B2 respectively, (b) The model reduces to a 
monolayer model with the third neighbor hopping t = t 2 /t± 
and the nearest neighbor hopping 2t (see the text). 



B2 sublattices to zero. The operators a\ - _(dj^ lCT ) rep- 
resent electron creation (annihilation) on the sublattice 
site Ai of the layer i with spin a =t,|, and b\ n a {bi 1 n,a) 
those for electrons on the sublattice site Bi. ijl is the 
chemical potential. We use units such that K=l. 



By 



introducing the Fourier transforms a r 



T,3 a i,l<r ex P(} k -3) and b i,k : a = £j^J, CT eX P( ife "J')' and 

diagonalizing the Hamiltonian, one obtains the spectrum 
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where a — 1 , 2 and ± denote 4 different branches of dis- 
persion and 



= i 



(3) 



In their orginal workpSl McCann and Fal'ko showed 
that the four-band model may be simplified to an ef- 
fective two-band model if one considers energies much 
smaller than t±. In momentum space, the Hamiltonian 
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in Eq. ([TJ) becomes 
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If we introduce the spinor 

^a(k) = (a 1>a ,k> a 2,<T,k> b 2,<r,k> b l,a,k) T ' 

the Hamiltonian can be expressed as a 4 x 4 matrix, 
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One may further define 2x2 matrices ffu = —fjj + 
t±a x , H 22 = -fil, #12 = -t(Re7 £ cr x +Im7gcr a ) = H 2 i, 
such that the eigenvalue equation can be written in the 
following form (fc indices are implied) 



#11 #12 
#21 #22 



*1 

* 2 



= E 



* 2 



from which we obtain 

{H 22 - ff 2 i (#n - E)- 1 H 12 }^2 = #*2- 



(8) 



(9) 



If we assume t± to be the largest energy scale (i.e. t± 3> /z 
) and apply the low energy limit (E ~ 0) , Eq. ([9]) becomes 
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*2 = #*2, 



(10) 



with tt 2 (A) = (6 2 ^,^^) T . 

As mentioned above, the two-band model described 
by the Hamiltonian in Eq. (10), is valid in the limit 
where E -C t± -C i. For energies larger than tj_j one 
needs to take into account the other two bands which 
overlap in energy with those considered in Eq. (10). 



If one takes into account trigonal warping at very low 
energies, the model is therefore valid in an energy win- 
dow lOmeV < E < 400 meV, where the dispersion is 
essentially isotropic and parabolic in the vicinity of the 
if -points, K± = 47rej / /3v / 3a. In the following sections, 
however, we are interested in superconducting order pa- 
rameters that account for the underlying lattice (point) 
symmetry beyond the parabolic approximation, and we 
use the simplified two-band model at even larger energies, 
up to the van-Hove singularity at ~ 3 eV. Formally, this 




FIG. 2. (Color online) Non-interacting dispersion (a) and 
density of states (b) of the projected monolayer model. Lin- 
ear dispersion in the vicinity of the K-points in the graphene 
monolayer (c) in comparison to the quadratic dispersion in 
our model (d). We use t = t 2 /t± for the unit of energy. 



amounts to increasing artificially the interlayer hopping 
t± such that it becomes the largest energy scale, t± 3> t. 
The Hamiltonian in Eq. ( 10 ) corresponds, in real 



space, to a single-layer honeycomb lattice with nearest- 
neighbor and third-neighbor hoppings. Whereas the ef- 
fective hopping amplitude of the latter is given by t 2 /t±, 
the effective nearest-neighbor hopping is twice as largePS 
This means that due to the strong interlayer hopping, 
the complete low-energy physics is projected onto the 
Bl and B2 sublattices which themselves form a hexago- 
nal lattice (see Fig. [l}. As mentioned above, the model 
is naturally valid in the small-momentum limit, i.e. for 
t 2 /t±\ka\ 2 ~ fj, -C t 2 /t± and reproduces correctly the fi- 
nite density of states (DOS) at E = of bilayer graphene 
[Fig.|. 

Since we consider the effective hopping t 2 /t± to be 
small and if there is a significant on-site repulsion U, 
spin-singlet bonds between Bl and B2 sites are expected 
to form due to superexchange processes. Therefore, we 
apply the t — J model but relax the requirement of the 
model that double occupation of sites is excluded. As we 
will be working in the mean-field approximation, we just 
assume an effective nearest-neighbor attractive interac- 
tion between electrons on Bl and B2 sublattices, and by 
doing this we favor spin-singlet bond formation. If this 
interaction is not too strong, it can be simply added to 
our Hamiltonian, with the help of the term 



H 



/ = -jVV& t - b,-> bK ^ , (ii) 



where J > 0. Now we apply the BCS ansatz by in- 
troducing the superconducting order parameter as a 3 
component complex vector 

A = (A^A^AsJ 
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FIG. 3. Different pairing instabilities in real space: (a) s- 
wave, (b) d x 2_ y 2 wave, (c) d xy wave, and (d) d x 2_ y 2 + id xy 
time reversal breaking d-wave. 



where the components are defined by 



(12) 



and correspond to the spin-singlet pairing amplitudes of 
three inequivalent pairs of nearest neighbors. The inter- 
action part Hi in the mean-field approximation becomes 



+27V^J|A S | 2 , (13) 
where TV is the number of unit cells. 



III. BOGOLIUBOV - DE GENNES ANALYSIS 
AND PAIRING SYMMETRIES 

The complete BCS Hamiltonian in momentum space 
is given by 



h = -— y (juk b,t +h.c 

k,a 
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H.c. 
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Similar to the case of the honeycomb monolayer,^ we 
can make our description much more transparent if we 



apply the following transformation that diagonalizes the 
kinetic part of the above Hamiltonian, 
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where = arg(7j). 

In this basis, where and represent the elec- 
tron states in the upper and lower band, respectively, 
the Hamiltonian transforms into 

H = 
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Here t = t 2 /t± and 
by 
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2 . The eigenvalues are given 



E- k = ±V (h) 2 + M 2 + 2J 2 (l^l 2 + IQI 2 ) ± 

(17) 

where C £ = A a cos(fc • u - 2<^), = £) a A 5 sin(fc • 
m — 2c^?g) and 

i4 = (^ 2 + 2J 2 |^| 2 )?e|+4J 4 (ReC*j:ImS , £ -ImCj:ReS' £ ) 2 . 

(18) 

If all Aa are purely real, i.e. there is no time-reversal 
symmetry breaking, then the second term in A is zero 
and the expression for the dispersion simplifies to 



2 + 2 J 2 52 
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(19) 



In this case S% only renormalizes the chemical poten- 
tial, whereas Cg plays the main role in the description of 
the superconducting order parameter. A comparison be- 
tween the Bogoliubov energy dispersion in Eq. ( 19 ) and 



the usual BCS expression shows that Cz can be identi- 
fied with the gap function. However, this name may be 
misleading because CV does not describe the gap, as in 
the example in Eq. (22 1 below. 



The symmetry analysis of the order parameter on a 
honeycomb lattice,^ yields the basis vectors which cor- 
respond to s, d x 2_ y 2 and d xy waves, respectively: 



1, 1) 



A(l, 

A= { A (2, -1, -1) 
A(0, 1, -1) 



(20) 



The three symmetries are illustrated in Fig. [3j and the 
function corresponding to these symmetries is shown 
in Fig. [4j in comparison with the monolayer case. The 
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FIG. 4. (Color online) Ck in the first Brillouin zone calculated 
for three possible symmetries on monolayer and projected bi- 
layer lattices. 



last two possibilities belong to a two-dimensional sub- 
space of irreducible representation of permutation group 
53!^ This means that any combination of these two order 
parameters is allowed by the symmetry. 

In the case of an s-wave order parameter with A = 
A (1,1,1), a small- wave- vector expansion (\q\a <C 1) 
around the X-points yields 



C K±+q 



V5 A 

T — q y aA. 



S K ± +q 



V3 



q x aA. (21) 



Thus both couplings are non-zero and no simple effective 
picture emerges by looking at the Hamiltonian in Eq.( 16 1. 



The lower excitation energy branch can be approximated 
in the small-momentum limit as 



(22) 



-[3 M t-(JA)2](|g»2, 



where we have used tf% . - 



9(k» 2 /4. 

If the coupling strengths are such that Eg has a mini- 
mum at q = 0, that is for (J A) 2 > 3 fit, a superconduct- 
ing instability may be realized. But as we pointed out 
in the introduction, an analysis based on the behavior in 
the region around q = is not well justified for the hon- 
eycomb bilayer if triagonal warping (linearly dispersing 
electrons) is not taken into account. If we nonetheless 
proceed with analysis we will get a time-reversal invari- 
ant superconducting instability with two kinds of Cooper 
pairs with p x +ip y and p x — ip y pairings. Due to the forms 

of Cr and Sr in the above Hamiltonian in the small mo- 
re k 

mentum limit, p-wave Cooper pairings are expected. For 
sufficiently large chemical potential [when we can neglect 
St in Eq.( 19 )], the system may be unstable towards a p y 



gapless superconductor, with gap minima on the Fermi 
surface, i.e. on a circle. 



For A = A(2, — 1, — 1), the small- momentum expan- 
sion around the X-points yields 



C K ±+ q{d x 2_ y 2) PS -3 
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and for A = A(0, 1, -1) 
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The gap function CV thus clearly shows the d x 2_ y 2 and 



the d xy symmetry in Eq. ( 23 ) and ( 24 ) , respectively. 



Notice that one may superpose two waves in the man- 
ner 



and 



C i {d±id) = C % {d x 2_ y 2) ± iV3C^(d xy ), (25) 



S%(d± id) = Sj:{d x 2^ y 2) ± iV3S^(d xy ), (26) 



which is identified with the d + id-wave superconducting 
phase in the following. In the small- wave- vector limit, 



the combined forms of C%, 
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restore the rotational symmetry - they are indeed eigen- 
states of rotation in two dimensions with the value of an- 
gular momentum equal to two. Only one form is allowed 
to exist in reality, and because it is the same irrespec- 
tive of the valley K or K' one obtains a solution that 
spontaneously breaks time-reversal symmetry. Thus we 
can identify the solution with the broken time-reversal 
symmetry d + id state. Something similar happens in 
the monolayer case, but the d-wave symmetry is recog- 
nized as a global dependence of the order parameter on 
the k vector in the Brillouin zone around central T-point 
(see Ref. I25p and p-wave behavior around K± points. 
In the bilayer case the time-reversal symmetry breaking 
d-wave order parameter emerges as a property of the low- 
energy small-momentum effective description around the 
K points, as shown above. 



IV. PHASE DIAGRAM 

We have found the ground state of our model Hamilto- 
nian for a broad range of J and p by minimizing the free 
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energy. At zero temperature, as a function of the order 
parameter, it is given by 

f=-EE %, Q + 2A^|A,| 2 , (29) 
feeiBZ a=±1 « 

where the first sum is over all wave vectors k in the 
first Brillouin zone and two Bogoliubov bands with pos- 
itive energies. The ground state is defined as a global 
minimum of the free energy in the order parameter 
space. In the present study, we concentrate on super- 
conducting order parameters in a variational approach, 
and thus we cannot exclude that other correlated (non- 
superconducting) phases may have an even lower energy. 
In the mean-field approach, superconducting ground 
states are expected even for infinitesimal positive values 
of J. 

The order parameter space is 6-dimensional, because 
it is defined by 3 complex numbers. However, adding 
the same phase to all three complex parameters does 
not modify the physical state, so one can always make 
one of the parameters purely real (we set a 2 real) 
and reduce the order parameter space dimensionality to 
5. We used the amoeba numerical method^ to directly 
minimize the free energy. Five-dimensional minimization 
often reveals more than one local minimum, but we were 
always able to identify the lowest-lying state to a sat- 
isfying level of certainty. However, for small values of 
J, the local free-energy minima are extremely shallow, 
with energies only slightly lower than the free energy of 
the normal state. Such features in the free-energy land- 
scape are completely clouded by numerical noise due to 
the discretization of the first Brillouin zone. Our numer- 
ical calculations are therefore limited to higher values of 
J, which give a solution with the amplitude of the order 
parameter larger than 10 -4 . This is marked by a dashes 
line on Fig. [5j and a short discussion of this line in the 
framework of a weak-coupling analysis may be found be- 
low [Eq. ((3TJ]. 

Our results are shown on Fig. [5] where the rele- 
vant quantities are represented by color in the (/x, J) 
plane. The amplitude of the order parameter is shown 
in Fig. 5(a). Upon small to moderate doping, the SC 
instability increases and becomes particularly favorable 
at the filling 5/8, which corresponds to the chemical po- 
tential fj,/t = 1, and the van-Hove singularity in the non- 
interacting DOS. For further doping the SC instability 
decreases. This gives to Fig. 5ja) roughly the look of the 
inverse DOS of Fig.gb). T fie gap in the single-particle 
excitations is shown in Fig. 5jh). It is particularly pro- 
nounced in the case of strong mixing of d x 2_ y 2 and id xy 
symmetry components, as we can see from Fig.[5|c). The 
contribution of different pairing symmetries is defined by 
the ratio w of different components of A, where 

A = A s e s + iA ls e s + A dxy e dxy + iA idxy e dxy 

+ A d x2 _ y2 e<i x2 _ y2 , (30) 

with e s = (1,1,1)^/3, e dxy = (0,1,-1),V2, and 
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FIG. 5. (Color online) (a) The order parameter amplitude, A, 
in the (/i, J) parameter space, obtained by a minimization of 
the free energy, (b) the single-particle excitation gap, (c) the 
contribution of id xy and (d) s-wave component in the ground 
state order parameter. The green dashed line marks where 
A drops below 10 -4 . Below this line, our numerics is not 
reliable. We use i = t 2 /t± for the unit of energy. 



e<l x 2_ 2 = (2, —1, — l)/v / 6- Fig-pic) shows the ratio 
w(id xy ) = \A idx J/\A\, and Fig. [Sp) the ratio w(s) = 
|A S |/|A|. The contributions of is and d xy components 
are negligible in all cases, and d x 2_ y 2 is the dominant 
component. 

The numerical results are, for clarity, also shown on 
Fig. [6] for three chosen values of the chemical potential, 
fj,/i = 0.04,0.55,1. Fig. [6|a) shows a sudden increase 
in the pairing amplitude with the increasing interaction 
J (note the logarithmic scale on the y-axis). For small 
J, the pairing amplitude is much larger for fi/t = 1, 
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FIG. 6. (Color online) (a) The order parameter amplitude A 
and (b) the single-particle excitation gap as a function of J, 
for fi = 0.04,0.55, 1. (c)-(e) The contributions of 3 relevant 
symmetry components. d x 2_ y 2 component is the dominant 
one for large J. The contribution of id xy increases with de- 
creasing J until the two contributions are equal and we find 
a pure d + id-wave symmetry. We use t = t 2 /t± for the unit 
of energy. 



i.e. at the van-Hove singularity, and in this case the 
single-particle excitation gap is also larger due the strong 
mixing of d x 2_ y 2 and id xy symmetries. Contributions 
of relevant components are compared in Figs. |6jc)-(e). 
At higher values of J one has a pure d x 2_ y 2 symmetry, 
whereas a mixture of d x 2_ y 2 and id xy symmetries is found 
at lower values of J. The contribution of id xy symmetry 
increases with decreasing J and almost pure d + id sym- 
metries are usually found at the lowest accessible values 
of J. 

Our numerical calculations were performed on proces- 
sors with 8GB of RAM which limited the number of k- 
points in the first Brillouin zone to 4000 x 4000, but we 
checked that results do not differ qualitatively even with 
a much sparser 2000 x 2000 fc-grid. A much denser and 
probably a non-uniform discretization of the first Bril- 
louin zone would be needed to probe the weak-coupling 
behavior of our model, that is for values of J below the 
dashed lines in Fig. [5] Notice, however, that the system 



in the small- J limit may be treated analytically within 
the weak-coupling limit the results of which are presented 
in Appendices A and B, for the cases of finite and zero 
chemical potential, respectively. 

In this weak-coupling regime and at finite chemical 
potential, we find that the d + id superconducting or- 
der parameter yields the lowest mean-field energy, when 
compared to order parameters that respect time-reversal 
symmetry (Appendix A), in agreement with our numer- 
ical results for larger values of J. Furthermore, this 
behavior is also expected from a theorem discussed in 
Ref. dnifor the BCS description of a single band, when 
both, d x 2_ y 2 and d xy , symmetries of the order parameter 
are allowed to exist. Although ours is a two-band (Bo- 
goliubov) model we expect a superconducting instability 
at any strength of attractive interaction at finite doping 
since the gap opens as 



J A ex exp 



8tt 1 



(31) 



(see Appendix A), in terms of the DOS p{Ep) at the 
Fermi level Ep. Notice, however, that the DOS is on 
the order of 1/t below the van-Hove singularity. There- 
fore, in order to resolve order parameters with a precision 
of 10~ 4 , as in our numerical calculations, the minimal 
coupling J ft should be on the order of unity, in rough 
agreement with the dashed line in Fig. [5j As already 
mentioned, the divergent DOS at the van-Hove singular- 
ity allows for a resolution of order parameters at even 
lower values of the coupling constant J. Equation (31 1, 



although derived in the weak-coupling limit, thus pro- 
vides a qualitative understanding of the above-mentioned 
inverse-DOS behavior of the order parameter in Fig. [5] 
Finally, we notice that the weak-coupling analysis 
yields a different picture at zero-doping (Appendix B), 
where a time-reversal-symmetric superconducting order 
parameter (with any real real combination of d x 2_ y 2 and 
d xy ) is energetically favored. 



V. RELEVANCE FOR BILAYER GRAPHENE 

In the following we discuss the relevance of our m odel 
for the physics of bilayer graphene. With an estimate^!! 
for the Coulomb on-site repulsion, U ~ 10 eV, intralayer 
nearest-neighbor hoppingpS t ~ 3 eV, and interlayer 
hoppingpSl £j_ ~ 0.4 eV, bilayer graphene has a tendency 
to develop strongly-correlated electron phases. Notice 
that, although similar energy scales are found in mono- 
layer graphene, the latter is to great accuracy described 
in terms of (quasi-)free electrons because of a vanish- 
ing DOS at the Fermi level, in the absence of intensive 
doping. 1 3 On the contrary, electronic correlations are 
much more efficient in bilayer graphene as a consequence 
of the finite DOS even at the band-contact points. This 
finite DOS may also be invoked when considering screen- 
ing. Whereas screening is highly inefficient in monolayer 
graphene, and one needs then to take into account the 



long-range nature of the electronic interaction potential, 
the screening properties in bilayer graphene are similar 
to those in usual 2D electron systems with a parabolic 
band dispersion. In this sense, an approach based on 
the Hubbard model, as used here excluding nearest and 
further-neighbor interactions, is better justified in bilayer 
than in monolayer graphene. Though in this sense ap- 
proximative our approach differs from other, mean-field 
and perturbative, approaches in the literature, which as- 
sume that interactions are weak despite their real mag- 
nitude. At half-filling we expect two antifcrromagnetic 

(Heisenberg) interactions, intra, J ~ jj ~ 1 eV, and in- 

t 2 

ter, Jj_ ~ ~ 16meV. Although the origin of the 
antiferromagnetic behavior may be claimed due to the 
quadratic dispersion of juxtaposed conduction and va- 
lence bancP^ (apart from triagonal warping at low en- 
ergies) the magnitude of U points out to also a role of 
Mott physics. Upon doping, holes or electrons will move 
on an effective (projected monolayer) lattice because in 
that way less antiferromagnetic bonds will be broken. 
Thus low energy physics will take place on the effective 
(B sites) lattice. An estimate from perturbative consid- 
erations for the effective antiferromagnetic coupling of B 
sites is J c ff ~ 2 = ~ 100 me V at finite t±. If we 

2 

rewrite the estimate for J c g as J c g ~ jj{jj) 2 = jjt\g, we 
can identify the effective hopping between Bl and B2 

.2 

sites to be t cS — ~ = J ~ l eV. Therefore the strongly- 
correlated regime of the bilayer graphene may be contin- 
uously related to our model by simultaneously increas- 

.2 

ing t± and decreasing U with t c g becoming t c g — ~, 

that is the hopping parameter t in our effective mono- 
layer model. If the on-site interaction is very large and 
precludes the double occupancy we can expect a situa- 
tion where doped holes or electrons are the only charge 

,2 

degrees of freedom and hop according t ~ f-j i.e. the 
hopping that characterize a non-interacting system (they 
do not pay the cost of U on their paths). We will as- 
sume that this is not the case in the graphene bilayer 
since U is not extremely large and keep our estimate of 

ieff = jj ~ 1 eV. 

Therefore modeled with two effective parameters, J c g 
and teg, bilayer graphene may be compared with the ef- 
fective honeycomb lattice considered in our paper and 
the corresponding t — J model. The main feature of bi- 
layer graphene appears to be that J c g ~ O.lteff <C t c g 
and in considering the relevance of our model we should 
confine ourselves to weak couplings, and small or moder- 
ate dopings; because we simplified the high-momentum 
physics of the bilayer (by considering the large t± limit) 
we should confine ourselves to lower dopings. First one 
sees from Fig. [5] that the gaps are in the meV range 
(2 to 5 meV for the maximal gaps) if one considers the 
energy scale t c g ~ J ~ I eV. Thus our results indicate 
very small energy scales that are unlikely to be resolved 
in today's graphene samples. Furthermore we should use 
t ff and J e ff for t and J for the exponent in the weak- 



coupling analysis in the Appendix A. Because we esti- 
mate t e g/ Jeff ~ 10, the weak-coupling analysis yields an 
exponential suppression and gaps below 1 meV, in agree- 
ment with our numerical findings shown in Fig. [5j 

VI. CONCLUSIONS 

In conclusion, we have investigated the possible oc- 
curence of d-wave superconductivity in doped bilayer 
graphene. We have used an effective lattice model in 
which the low-energy electrons are allowed to hop on a 
honeycomb lattice that consists of the B sites of the two 
layers. This model reproduces the low-energy physical 
properties of bilayer graphene below t± ~ 400 meV and 
neglects the effect of higher-energy bands above this value 
by artificially increasing the value of t± , while maintain- 
ing the lattice point symmetry. This point symmetry 
happens to be important in fixing the symmetry of the su- 
perconducting order parameter when increasing the cou- 
pling constant J > for a bond-singlet order within a 
t — J model that has been investigated in a mean-field 
approximation. 

At low values of J, a time-reversal symmetry-breaking 
d + id superconducting order may be sta bilized, s imi- 
larly to the propo sed phases in moderately^EsHSl anc j 
highly dopecP^ monolayer graphene, such as to re- 
store rotational symmetry. The time-reversal symme- 
try breaking in this phase is apparent in the effective 
(low-momentum) description ~ (k x + ik y ) 2 at both val- 
ley points, K and K' . Most saliently, the low- J d + id- 
wave superconducting order persists upon doping until 
the van-Hove singularity at 5/8 filling for electron dop- 
ing (or 3/8 for hole doping). For larger values of J, 
rotational symmetry is no longer respected because of 
the increased importance of lattice point symmetry. The 
anisotropic time-reversal symmetric d x 2_ y 2-wave compo- 
nent then becomes increasingly important in the order 
parameter, such that more conventional d-wave super- 
conductivity may be expected. One notices, however, 
that the superconducting energy scales, estimated both 
numerically and within a weak-coupling analysis, are very 
small (in the meV range), such that an experimental ob- 
servation in presently available graphene samples is un- 
likely, and an important increase in sample quality is 
required. Similarly to monolayer graphenej^ one may 
though expect a renormalization flow to stronger cou- 
pling constants when approaching the van-Hove singu- 
larity in bilayer graphene by chemical doping. 
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Appendix A: Weak-coupling analytical solution at 
finite chemical potential 

Here, we present briefly the weak-coupling analysis of 
superconducting order in the effective bilayer model. In 
order to simplify the notation, we use the letter t to de- 
note the effective hopping t. The DOS at the Fermi level, 
p{Ep), is on the order of the inverse hopping parameter 
l/t. Notice that, if only a parabolic band is taken into 
account it remains fixed at its Ep = value, but correc- 
tions to the parabolic approximation immediately yield 
a contribution that varies linearly with the Fermi level, 
in agreement with the DOS plotted in Fig. [2jb). 

In the case when A = A(l, 1, 1), a weak-coupling BCS 
analysis that takes into account only electrons in the 
lower Bogoliubov band gives 



Because of its twice lower mean-field energy, the d x 2_ y 2 + 
i\/3d xy time-reversal symmetry breaking instability, 
which we call in short d-wave, is more likely than d x 2_ y 2 
and d xy -w&ve order parameters. In the large-doping 
limit, the energy minimization is also much more effi- 
cient for d-wave than p y - wave as seen in the small value 
of the ratio 



6E P MF 

6E tlF 



JJ_ 

2E r 



exp 



2?r x 8 



P(E F )J 



9t 

2fJL 



The most natural choice for E c 



1 



(A7) 
is to be of 



for ji < y 

the order of as a first energy scale when we start from 
the smallest one, i.e. J. The time-reversal symmetry 
breaking d-wave is the solution of our BCS mean-field 
Hamiltonian also due to the theorem proved in Ref. I2U1 
In the case of weak coupling that we consider here, i.e. 
J <C /i, and p > (electron doping), we have an effective 
one-band theory of upper-band electrons to which the 
theorem can be applied. 

In the following we will look more closely into an ef- 
fective, low-energy description of the d- wave instability 
in the case of high electron doping and discuss only the 
lower energy Bogoliubov band. Therefore our effective 
Hamiltonian is 



J A = y/2tE c exp -24V3tt 



pp(E F )J 



(Al) 



with E c as an energy cut-off around the Fermi value, for 
the solution, and 



SEl 
N 



wf = _ (JA) 2W(M 1 , (A2) 



t 4V3V 



for the gain in the mean-field energy, 5Emf, by the pair- 
ing instability. 

The weak coupling BCS analysis in the case of electron 
doping (p > 0) for d x 2_ y 2 and d x 2_ y 2 + i\/3 d xy gives 

for the solution which we denoted by A = A^, and 

' a =#M-;w+0' ,A4) 

in the case of d xy wave. For the energy gain one obtains 



5E M F{d x 2- y 2) _ SE MF (d. 



■r:yj 



N 



N 



-(JA d ) 2 p(E F ) 



and for a d, r 2_„2 + iyo d xy wave, one finds 



3V3 
(A5) 



SE 



MF -^-(JA d ) 2 p(E F )^. (A6) 



He = B%-^) c l c ^+£ (VIA + H - c -) (A8) 

kcr k 



where A^ 



N 



2tt 



(k x — iky) 2 /\k\ 2 . In the weak-coupling BCS 
analysis it can be easily shown that the Hamiltonian is 
completely equivalent to the one with A^ ~ (k x — ik y ) 2 , 
because both Hamiltonians have an effective description 
on a Fermi circle defined by ie^ = p. With this ad- 
justment we have exactly the form of the BCS Hamil- 
tonian studied in Ref. [35] on time-reversal symmetry 
breaking superconductors in two dimensions. In the so- 
called weak-pairing case for finite p > that we want to 
study, the minimum of Bogoliubov excitations moves to 
finite values of k, te^ = p, i.e. to the Fermi surface of 
free particles. The Cooper pair wave function g(r) may 
be a non-universal function of \r\ where r is the relative 
coordinate of the pair. On the other hand, the depen- 
dence of the function on the angle of vector r is fixed 
and can easily be derived in the Bogoliubov formalism 
to be 3(H) oc | cx (x — iy) 2 where z — x + iy is the 
two-dimensional complex coordinate. Thus the relative 
angular momentum of the Cooper pair is I = —2. The 
weak-pairing phase is topological, gapped in the bulk be- 
cause p > 0, and possesses a doublet of spin 1/2 Dirac 
edge modes^S. In our case, because of the fermion dou- 
bling on the honeycomb lattice and the existence of the 
two K points (valleys) [and because around each one we 
have the same effective description given by Hamiltonian 
in Eq.(16l], we expect four Dirac modes on the edge. 
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Appendix B: Weak coupling analytical solution at 
zero chemical potential 



for the energy gain. On the other hand for d + id sym- 
metry we find 



In the weak coupling limit at ji = 0, when both Bo- 
goliubov bands are taken into account we find for d x 2_ y 2 
symmetry 



JA 



-5c 



2c 



and 



(B3) 



JA d = ^ exp 



He 

2c 



fBll 



Because 



= -9c (JA d+id ) 2 . 



(B4) 



with c : 



271-^3 * 



i for the solution, and 



r rpd+id 



M F UJ2j MF 



S rpd+id 



(B5) 



di ^ = -^c(JA rf ) 2 , 



any real combination of d x i_ y i and d^j, waves is more 
(B2) likely than d + id wave. 
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